; docformat = 'rst'
;
;+
;
; :Purpose:
;   Define LPDM and EUM domains and grids, based on values given in 
;   {input_directory}/{sim_name}/EUM/grid.inp.sav and Lon_domain and Lat_domain parameters.
;   
; :Inputs:
;   sim_name: Simulation name
;
; :Requires:
;   {input_directory}/{sim_name}/EUM/grid.inp.sav file
;   
;   LON_DOMAIN input parameter ({sim_name}.params file)
;   
;   LAT_DOMAIN input parameter ({sim_name}.params file)
;   
;   STATIONS input parameter ({sim_name}.params file)
;
; :Outputs:
;   {input_directory}/regions.sav
;
; :History:
; 	Written by: Matt Rigby, MIT, Aug 19, 2011
;
;-
pro cels_define_regions, sim_name

  compile_opt idl2
  on_error, 2

  lon_domain=cels_get_parameter(sim_name, 'LON_DOMAIN')
  lat_domain=cels_get_parameter(sim_name, 'LAT_DOMAIN')
  stations=cels_get_parameter(sim_name, 'STATIONS')

  N_LR=n_elements(lon_domain[0, *])

  ;Get EUM grid
  restore, cels_filestr(/Input, sim_name + '/EUM/grid.inp.sav')
  Lon_EUM=lon
  Lat_EUM=lat

  lon_domain_EUM=lon_domain
  lat_domain_EUM=lat_domain

  dlon=lon_EUM[1] - lon_EUM[0]
  dlat=lat_EUM[1] - lat_EUM[0]

  for si=0, N_LR-1 do begin
    lon_EUM_temp=lon_EUM
    if lon_domain[0, si] lt 0. then begin
      lon_EUM_temp[where(lon_EUM gt 180.)]=lon_EUM[where(lon_EUM gt 180.)] - 360.
    endif

    dummy=min(abs(lon_EUM_temp - lon_domain[0, si]), LonM)
    dummy=min(abs(lon_EUM_temp - lon_domain[1, si]), LonP)   
    lon_domain_EUM[0, si]=lon_EUM[lonM]
    lon_domain_EUM[1, si]=lon_EUM[lonP]
    lon_domain[0, si]=lon_EUM_temp[lonM]
    lon_domain[1, si]=lon_EUM_temp[lonP]

    dummy=min(abs(lat_EUM - lat_domain[0, si]), LonM)
    dummy=min(abs(lat_EUM - lat_domain[1, si]), LonP)   
    lat_domain_EUM[0, si]=lat_EUM[lonM]
    lat_domain_EUM[1, si]=lat_EUM[lonP]
  endfor
  lat_domain=lat_domain_EUM
  

  nStations=n_elements(stations)
  Station_LR=intarr(nStations)

  restore, cels_filestr(/Input, sim_name + '/measurements/measurements.sav')

  ;Find which stations are in which local region (defined by lon_domain and lat_domain)
  for si=0, nStations-1 do begin
    wh_station=where( mf_station_ID eq stations[si], count )
    if count ne 1 then begin
      message, 'Station missing in measurement.sav, or defined multiple times', stations[si]
    endif
    station_LR[si]=mr_grid_inside(mf_lon[wh_station], mf_lat[wh_station], lon_domain, lat_domain)
  endfor
  
  ;Define local grids based on footprint files
  ;Assumes that if there are multiple stations are in the same region, they use the same grid
  ;Footprints are assumed to be stored in annual files
  lon_local=hash()
  lat_local=hash()
  
  For ri=0, N_LR-1 do begin
    wh=where(station_LR eq ri, count)
    if count eq 0 then message, 'Broken' ;SOMETHING STRANGE HAS HAPPENED

    fname=!null
    for whi=0, n_elements(wh)-1 do begin
      fname=[fname, file_search(cels_filestr(/Input, 'fp/'), strcompress('*' + stations[wh[whi]] + '*', /remove_all))]    
    endfor

    if fname eq !null then begin
      message, "Can't find footprint files"
    endif

    restore, fname[0]
    lon_local[ri]=lon
    lat_local[ri]=lat
    ;Check that other stations in region have same grid
    if n_elements(fname) gt 1 then begin
      for fi=0, n_elements(fname)-1 do begin
        restore, fname[fi]
        if (n_elements(lon_local[ri]) ne n_elements(lon)) or (n_elements(lat_local[ri]) ne n_elements(lat)) then begin
          message, strcompress("Footprint grids different dimension: " + fname[fi])
        endif
        if (total(abs(lon_local[ri] - lon)) gt n_elements(lon)*(lon[1] - lon[0])/10.) or $
          (total(abs(lat_local[ri] - lat)) gt n_elements(lat)*(lat[1] - lat[0])/10.) then begin
          print, strcompress("CELS_DEFINE_REGIONS: WARNING, footprint grids different locations: " + fname[fi])
        endif
      endfor
    Endif    
  endfor


  save, filename=cels_filestr(/Input, sim_name + '/regions.sav'), lon_domain, lat_domain, lon_domain_EUM, lat_domain_EUM, station_LR, $
    N_LR, lon_EUM, lat_EUM, lon_local, lat_local

end